﻿import numpy as np
import geatpy as ea
import problem1
import pandas as pd
import networkx as nx
from matplotlib.font_manager import FontProperties

# 加载中文字体
font = FontProperties(fname="./Deng.ttf", size=14)


class MyProblem(ea.Problem):  # 继承Problem父类
    def __init__(self):
        name = 'MyProblem'  # 初始化name（函数名称，可以随意设置）
        M = 1  # 初始化M（目标维数）
        maxormins = [1]  # 初始化maxormins（目标最小最大化标记列表，1：最小化该目标；-1：最大化该目标）
        self.var_set = np.array([1.1, 1, 0, 3, 5.5, 7.2, 9])  # 设定一个集合，要求决策变量的值取自于该集合

        # 共120个决策变量，前24个变量为24个车的出发点等待时间，后面的变量依次为24个第一次装载点、24个第一次发射点、24个第二次装载点、24个第二次发射点
        Dim = 120  # 初始化Dim（决策变量维数）

        # 前24个变量为24个车的出发点等待时间，可以是小数。后面都是用整数表示各个地点，地点和整数对应关系如下
        # 地点 D1,D2, Z01,Z02,...,Z06, F01,F02,...,F60, J01,J02,...,J62
        # 整数 0, 1 , 2  ,3  ,...,7  , 8  ,9  ,...,67 , 68 , 69,...,129
        # 另外，每个车也对应一个整数
        # 车   A1,A2,A3,A4,A5,A6, B1,...,B6, C1,...,C12
        # 整数  0,1 ,2 ,3 ,4 ,5 , 6 ,...,11, 12,...,23
        varTypes = [0] * 24 + [1] * (Dim - 24)  # 初始化varTypes（决策变量的类型，元素为0表示对应的变量是连续的；1表示是离散的）
        self.z_points_start, self.z_points_end = 2, 7  # 见上表就知道Z01对应整数为2，Z06对应整数为7
        self.f_points_start, self.f_points_end = 8, 67  # 见上表就知道F01对应整数为8，F60对应整数为67
        # 至于J01～J62，这里我们先不加入到决策变量中，因为J只是路过的节点而已，不是目标函数需要关心的
        # 我们通过遗传算法计算就能得到经过的J，也就是每个车的路径了

        # 下面这些是决策变量的上下界
        t_lb, t_up = [0.0] * 24, [60.0] * 24  # 人为设定出发点等待时间在0.0小时到60.0小时之间
        z_lb, z_up = [self.z_points_start] * 24, [self.z_points_end] * 24  # 24辆车对应的装载点对应整数范围
        f_lb, f_up = [self.f_points_start] * 24, [self.f_points_end] * 24  # 24辆车对应的发射点对应整数范围
        lb = t_lb + z_lb + f_lb + z_lb + f_lb  # 决策变量下界
        ub = t_up + z_up + f_up + z_up + f_up  # 决策变量上界
        lbin = [1] * Dim  # 决策变量下边界（0表示不包含该变量的下边界，1表示包含）
        ubin = [1] * Dim  # 决策变量上边界（0表示不包含该变量的上边界，1表示包含）
        # 调用父类构造方法完成实例化
        ea.Problem.__init__(self, name, M, maxormins, Dim, varTypes, lb, ub, lbin, ubin)

        # 读点数据进来，用于得到邻接矩阵
        pos = pd.read_excel("data.xls", sheet_name=0, index_col=0)
        connections = pd.read_excel("data.xls", sheet_name=1, header=None)

        # 下面这部分求解图论的内容，具体为得出邻接矩阵，解出每个地点两两的最近路径矩阵
        # self.neighbour就是邻接矩阵，就是130行130列，下面列出了行和列对应的地点，x就是实际数据，...就是写不下省略
        #     D1,D2, Z01,Z02,...,Z06, F01,F02,...,F60, J01,J02,...,J62
        #  D1 x  x   x   x   ...  ... ... ... ... ... ... ... ...  x
        #  D2 x  x   x   x   ...  ... ... ... ... ... ... ... ...  x
        # Z01 x  x   x   x   ...  ... ... ... ... ... ... ... ...  x
        # Z02 x  x   x   x   ...  ... ... ... ... ... ... ... ...  x
        # ... ... ... ... ... ... ... ... ... ... ... ... ... ... ...
        # Z06 ... ... ... ... ... ... ... ... ... ... ... ... ... ...
        # F01 ... ... ... ... ... ... ... ... ... ... ... ... ... ...
        # ... ... ... ... ... ... ... ... ... ... ... ... ... ... ...
        # F60 ... ... ... ... ... ... ... ... ... ... ... ... ... ...
        # J01 ... ... ... ... ... ... ... ... ... ... ... ... ... ...
        # ... ... ... ... ... ... ... ... ... ... ... ... ... ... ...
        # J62 x  x   x   x   ...  ... ... ... ... ... ... ... ...  x
        self.neighbour = problem1.convert_neighbour(connections, pos)
        graph = nx.from_numpy_matrix(self.neighbour)  # 用邻接矩阵建立无向有环图
        # self.shortest_paths指每个地点两两的最近路径矩阵
        # 和邻接矩阵一样的排列方式，只不过每个数据的位置是一个列表，这个列表表示从所在行对应地点到所在列对应地点的最短路径
        # 比如Z01（对应整数2）到Z06（对应整数7）最短路径列表是self.shortest_paths[2][7]=[Z01,Z02,Z03,Z04,Z05,Z06]
        self.shortest_paths = [problem1.shortest_path(i, graph) for i in range(130)]

        # 下面这部分指定24个车各自的的起始点和各自的速度，也指定了路段的类型，速度分两个，单向路速度和双向路速度
        # 通过公式来确定某条路上某个车的速度，如下
        # speed = (1-self.road_type[路段端点1地点对应的整数, 路段端点2地点对应的整数]) * self.car_speed_n[车对应的整数] +
        #          self.car_speed_d[车对应的整数]*self.road_type[路段端点1地点对应的整数, 路段端点2地点对应的整数]
        # 车均匀分布在两起始点  A1 A2 A3 A4 A5 A6 B1 B2 B3 B4 B5 B6 C1 C2 C3 C4 C5 C6 C7 C8 C9 C10C11C12
        self.start_point = [0, 0, 0, 1, 1, 1, 0, 0, 0, 1, 1, 1, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1]  # 0是D1对应整数，1是D2对应整数
        #  单向道路的车速      A1  A2  A3  A4  A5  A6  B1  B2  B3  B4  B5  B6  C1  C2  C3  C4  C5  C6  C7  C8  C9  C10 C11
        #  后缀为n(normal)   C12
        self.car_speed_n = [45, 45, 45, 45, 45, 45, 35, 35, 35, 35, 35, 35, 30, 30, 30, 30, 30, 30, 30, 30, 30, 30, 30,
                            30]
        #  双向道路的车速      A1  A2  A3  A4  A5  A6  B1  B2  B3  B4  B5  B6  C1  C2  C3  C4  C5  C6  C7  C8  C9  C10 C11
        #  后缀为d(double)   C12
        self.car_speed_d = [70, 70, 70, 70, 70, 70, 60, 60, 60, 60, 60, 60, 50, 50, 50, 50, 50, 50, 50, 50, 50, 50, 50,
                            50]
        # 道路类型 self.road_type与邻接矩阵形状形式一致，其中数据0代表单向道路，1代表双向道路, inf代表没有直接连接
        self.road_type = self.neighbour.copy()
        self.road_type[self.neighbour < np.inf] = 0  # 邻接矩阵中有直接连接的先初始化为单行道
        self.road_type[68:87, 68:87] += 1  # J01（对应整数68）到J11以及J12到J20（对应整数87）为双向的
        self.road_type = np.asarray(self.road_type, dtype=np.uint8)

        # 用于记录每个点的使用次数，只不过存储的是整数

    def aimFunc(self, pop):  # 目标函数
        Vars = pop.Phen  # 行数为个体数，列数就是120个决策变量，前24个变量为24个车的出发点等待时间，后面的变量依次为24个第一次装载点、24个第一次发射点、24个第二次装载点、24个第二次发射点
        obj_v = np.zeros((Vars.shape[0], 1))
        cv = np.zeros((Vars.shape[0], 1))
        for geti_idx in range(Vars.shape[0]):
            # 130个元素的向量，分别代表D1,D2, Z01,Z02,...,Z06, F01,F02,...,F60, J01,J02,...,J62的使用次数
            position_repeat = np.zeros((130,))
            # solution整成5行24列，如下列数为车数，行数为出发等待时间和各地点
            #                    A1,A2,A3,A4,A5,A6, B1,...,B6, C1,...,C12
            # 出发点等待时间/小时    x  x  x  x  ...  ... ... ... ... ... x   <-- solution_time_wait[:,:]
            # 第一次装载点对应整数   x  x  x  x  ...  ... ... ... ... ... x    <-- solution_route[0,:]
            # 第一次发射点对应整数   x  x  x  x  ...  ... ... ... ... ... x    <-- solution_route[1,:]
            # 第二次装载点对应整数   x  x  x  x  ...  ... ... ... ... ... x    <-- solution_route[2,:]
            # 第二次发射点对应整数   x  x  x  x  ...  ... ... ... ... ... x    <-- solution_route[3,:]
            # 只看首行，就是solution_time_wait。去除首行，就是solution_route
            solution = Vars[geti_idx, :].reshape((5, 24))
            solution_time_wait = np.asarray(solution[0, :], dtype=np.int)
            solution_route = np.asarray(solution[1:, :], dtype=np.int)
            # routine用来形成如下的列表，一行就是一个车的两次发射路线，因此每一行列数不一定相等
            #     出发点 中间地点 中间地点 ...  第一次装载点 中间地点 ... 第一次发射点 中间地点 ... 第二次装载点 中间地点 ... 第二次发射点
            # A1  0     x      x       ...  x         x       ... x         x      ...  x         x       ... x
            # A2  0     x      x       ...  x         x       ... x         x      ...  x         x       ... x
            # ... ...   ...    ...     ...  ...       ...     ... ...       ...    ...  ...       ...     ... ...
            # A6  1     x      x       ...  x         x       ... x         x      ...  x         x       ... x
            # B1  0     x      x       ...  x         x       ... x         x      ...  x         x       ... x
            # ... ...   ...    ...     ...  ...       ...     ... ...       ...    ...  ...       ...     ... ...
            # B6  1     x      x       ...  x         x       ... x         x      ...  x         x       ... x
            # C1  0     x      x       ...  x         x       ... x         x      ...  x         x       ... x
            # ... ...   ...    ...     ...  ...       ...     ... ...       ...    ...  ...       ...     ... ...
            # C12 1     x      x       ...  x         x       ... x         x      ...  x         x       ... x
            routine = []
            # time最终这样形状，表示各车各路段所用时间，小时为单位
            #                    A1,A2,A3,A4,A5,A6, B1,...,B6, C1,...,C12
            # 出发点->中间地点      x  x  x  x  ...  ... ... ... ... ... x
            # 中间地点->...        x  x  x  x  ...  ... ... ... ... ... x
            # 第一次装载点->中间地点 x  x  x  x  ...  ... ... ... ... ... x
            # ...                ... ... ... ... ... ... ... ... ... ...
            # 最后第二行举例        0   0   0   0   0   x   x   0   0   x
            # 最后一行举例         0   0   0   0   0   0   x   0   0   x
            # 由于不知道各车会经过几个路段，所以行数为最长的那个车经过的路段数。其余车不足的地方补0
            # 不用担心这个999，后续会截断到正确的行
            time = np.zeros((999, 24))
            # i表示车对应的整数
            for i in range(24):
                start_point = self.start_point[i]  # 取得当前车（对应整数i）的出发点对应的整数
                routine.append(self.shortest_paths[start_point][solution_route[0, i]])  # 获取当前车辆（对应整数i）从出发点到第一次装载点的最短路径
                for j in routine[-1][:-1]:
                    position_repeat[j] = position_repeat[j] + 1
                    pass
                node_idx_for_time = 0
                # j表示个der（第一声）
                for j in range(3):
                    # 在循环里依次获取获取当前车辆的r：
                    # 从第一次装载点到第一次发射点的最短路径   <-- j=0
                    # 从第一次发射点到第二次装载点的最短路径   <-- j=1
                    # 从第二次装载点到第二次发射点的最短路径   <-- j=2
                    r = self.shortest_paths[solution_route[j, i]][solution_route[j + 1, i]]
                    routine.append(r)  # 将r记录到routine里

                    # 遍历r
                    for k in range(len(r) - 1):
                        p1, p2 = r[k], r[k + 1]
                        # 算速度，套了前面注释的的公式
                        speed = (1 - self.road_type[p1, p2]) * self.car_speed_n[i] + self.car_speed_d[i] * \
                                self.road_type[
                                    p1, p2]
                        # 算时间，单位小时，记录到time里
                        time[node_idx_for_time, i] = self.neighbour[p1, p2] / speed
                        node_idx_for_time = node_idx_for_time + 1

                        position_repeat[p1] = position_repeat[p1] + 1
                position_repeat[routine[-1][-1]] = position_repeat[routine[-1][-1]] + 1

            # time矩阵为24列n行的矩阵。24列表示有24辆车
            # 比如A2车从D1（序号为0）出发，到J10（77），到J09（76），到J45（112）...
            # 那么time[:,1]=[[xxx], A2车从D1到J10路段用的时长
            #               [xxx], A2车从J10到J09路段用的时长
            #               [xxx], A2车从J09到J45路段用的时长
            #               ...]
            time = time[~(time == 0).all(1), :]

            # 对上面的各路段时长做累加操作，并认为在0时刻出发，例子同上
            # time_acc[:,1]=[[xxx], A2车从D1出发到J10的时刻
            #               [xxx], A2车到J09的时刻
            #               [xxx], A2车到J45的时刻
            #               ...]
            time_acc = np.add.accumulate(time, 0)  # 首行各车的出发等待时间这里没有加上

            # 当然出发前其实有等待，不能认为0时刻出发，因此加上等待时间，注意这里的等待时间
            # time_acc_with_wait = time_acc + solution_time_wait

            # time_acc最后一行就是所有车到终点各自用掉的时间，加起来就是暴露时间#######################此处存疑 待确认#########################
            obj_v[geti_idx, 0] = np.sum(time_acc[-1, :])
            cv[geti_idx, 0] = np.max(position_repeat[self.f_points_start:self.f_points_end + 1]) - 1
        pop.ObjV = obj_v
        pop.CV = cv


if __name__ == '__main__':
    """===============================实例化问题对象==========================="""
    problem = MyProblem()  # 生成问题对象
    """=================================种群设置=============================="""
    Encoding = 'RI'  # 编码方式
    NIND = 20  # 种群规模
    Field = ea.crtfld(Encoding, problem.varTypes, problem.ranges, problem.borders)  # 创建区域描述器
    population = ea.Population(Encoding, Field, NIND)  # 实例化种群对象（此时种群还没被初始化，仅仅是完成种群对象的实例化）
    """================================算法参数设置============================="""
    myAlgorithm = ea.soea_DE_rand_1_bin_templet(problem, population)  # 实例化一个算法模板对象
    myAlgorithm.MAXGEN = 25  # 最大进化代数
    myAlgorithm.mutOper.F = 0.5  # 差分进化中的参数F
    myAlgorithm.recOper.XOVR = 0.2  # 重组概率
    myAlgorithm.trappedValue = 1e-6  # “进化停滞”判断阈值
    myAlgorithm.maxTrappedCount = 10  # 进化停滞计数器最大上限值，如果连续maxTrappedCount代被判定进化陷入停滞，则终止进化
    myAlgorithm.logTras = 1  # 设置每隔多少代记录日志，若设置成0则表示不记录日志
    myAlgorithm.verbose = True  # 设置是否打印输出日志信息
    myAlgorithm.drawing = 1  # 设置绘图方式（0：不绘图；1：绘制结果图；2：绘制目标空间过程动画；3：绘制决策空间过程动画）
    """===========================调用算法模板进行种群进化======================="""
    [BestIndi, population] = myAlgorithm.run()  # 执行算法模板，得到最优个体以及最后一代种群
    BestIndi.save()  # 把最优个体的信息保存到文件中
    """==================================输出结果=============================="""
    print('用时：%f 秒' % myAlgorithm.passTime)
    print('评价次数：%d 次' % myAlgorithm.evalsNum)
    if BestIndi.sizes != 0:
        print('最优的目标函数值为：%s' % BestIndi.ObjV[0][0])
        print('最优的控制变量值为：')
        for i in range(BestIndi.Phen.shape[1]):
            print(problem.var_set[BestIndi.Phen[0, i].astype(int)])
    else:
        print('没找到可行解。')
